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We investigate dispersive optical nonlinearities that arise from Rydberg excitation blockade in 
cold Rydberg gases. We consider a two-photon transition scheme and study the non- linear response 
to a weak optical probe in presence of a strong control beam. For very low probe fields, the 
dominant nonlinearities are of the third order and they can be exactly evaluated in a steady state 
regime. In a more general case, the change in average atomic populations and coherences due 
to Rydberg interactions can be characterized by properly defined scaling parameters, which are 
generally complex numbers but in certain situations take the usual meaning of the number of atoms 
in a blockade sphere. They can be used in a simple "universal scaling" formula to determine the 
dispersive optical nonlinearity of the medium. We also develop a novel technique to account for the 
Rydberg interaction effects, by simplifying the treatment of nonlocal interaction terms, the so-called 
collisional integrals. We find algebraic relations that only involve two-body correlations, which can 
be solved numerically. All average populations and coherences are then obtained straightforwardly. 



I. INTRODUCTION 

The strong interactions between Rydberg atoms pQ 
are studied and exploited in various fields ranging from 
molecular and ultracold plasma physics [2j |3] to nonlin- 
ear optics and quantum computing [U [5]. In particu- 
lar, many efforts are currently dedicated to create large 
optical non-linearities in Rydberg atomic ensembles and 
make them observable at a few-photon level in order to 
use them for quantum information processing. 

Using an atomic cloud as a non-linear optical medium 
often involves Electromagnetically Induced Transparency 
(EIT) schemes [6 , where the propagation of one beam 
(probe) is strongly modified by the presence of another 
beam (control), due to a two-photon resonance condition. 
By varying the control field, one can convert optical ex- 
citations into atomic ones, creating dark-state polaritons 
[7] which propagate through the atomic cloud at a very 
slow speed. This increases the non-linear effects related 
to the anharmonic level structure of each individual atom 
[8HIU]. but so far it was not sufficient to make them ob- 
servable at a few photon level. On the other hand, it has 
been demonstrated that when the two-photon transition 
involves a Rydberg state, the non-linear response of the 
atomic cloud can be strongly enhanced [4j [TTHT4] . In a 
simple physical picture, this happens because the inter- 
actions bring neighboring atoms out of the two-photon 
resonance, so each excited Rydberg atom becomes sur- 
rounded by a blockade volume where no other Rydberg 
excitation is possible. These "blockade" interactions be- 
tween the polaritons increase the optical non-linearity 
which can become significant even at a single-photon level 
[15j[T6]. In the resonant, dissipative regime, it has re- 
cently been demonstrated that this phenomenon could 
be used to operate the atomic cloud as a single-photon 
source [T7U22]. 

Many applications in quantum information process- 
ing require to transpose these effects from the dissipa- 
tive to the dispersive regime, in order to enable coherent 
photon-photon interactions. Large Rydberg-induced dis- 



persive nonlinearities have been recently measured [13] 
even though the regime of single-photon nonlinearities 
has not been reached. In this paper we study the dis- 
persive nonlinearities in a Rydberg EIT medium using 
several theoretical methods, with experimental parame- 
ters corresponding to ref. [13 . The results are compared 
to the experiment, and the advantages and disadvantages 
of these various methods are discussed in detail. 



II. THEORY 

We consider an ensemble of TV ^> 1 three-level atoms 
with a ground state |1), a short-lived intermediate state 
1 2), and a highly excited long-lived Rydberg state |3) 
(see Fig. [I]). We study the non-linear response of the 
ensemble to a probe beam with a Rabi frequency Q p and 
optical frequency u p , detuned from the |1) — >> |2) transi- 
tion by A2 = uo p — CJ12, in presence of a strong control 
beam with a Rabi frequency Q c and optical frequency co> c , 
with a two-photon detuning A3 = uj p + uo c — 6J13. It is 
assumed that the interaction between two atoms i and j 
in a Rydberg state is described by a single potential kij 
(as discussed in Section 



IV 



in case of multiple potential 
curves one may use an effective single potential). In the 
rotating-wave approximation, the many-body dynamics 
is governed by the Hamiltonian 
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where & % a ^ = |c^)(/^|, and Q c is taken real. We assume 
that the Rydberg ensembles are locally uniform meaning 
that the spatial variations of any a l a ^ and of the atomic 
density rj are on a significantly larger scale than the cor- 
relation length imposed by the interactions. Because of 
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FIG. 1. a) An ensemble of three- level atoms is excited by a 
strong (blue) control field, and a weak (red) probe field. The 
two-photon detuning A3 is much less in magnitude than the 
single photon detuning A2. b) Effects of Rydberg interactions 
in the vicinity of the two-photon resonance. Since the strong 
interactions effectively decouple the atoms from the control 
field, the blockaded atoms behave as two-level atoms which 
means that the resonance feature (blue) of the susceptibility 
of three-level non-interacting atoms gradually diminishes and 
moves towards the featureless two-level dependence. 



this, we can drop all the indices z, j unless keeping them 
is necessary for the clarity of the expressions. 

To include the decay of the intermediate state and fi- 
nite laser linewidths, relaxation constants T a p = j a p — 
i(Ap — A a ) associated to each a a /3 are introduced. Here 
7 a/ 3 are the total (positive) decay constants of the co- 
herences, and Ai = 0. The optical Bloch equations for 
single-atom averages a l a a — are 
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where V£p = ^j^i^ij^a^ 33) are two-atom collisional 
integrals. In the following, we use the shorter notation 
= (° % olP°Ilv) and assume the symmetry relation 
Wa(3,ij,v — ov^ap. For numerical evaluations, we will 
express all frequencies in units of the linewidth 7 = 712 = 
722/2 of the |1) — >• 1 2) transition. 

For some general considerations, it may be useful to 
solve Eqs. for the averages a a p assuming that the 

collisional integrals are known. Obviously, in the steady 
state such solutions are linear combinations of the colli- 
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sional integrals V 13 , V 3 i, V 23 , and V 32 
^ = tf™ + P%V 13 + 0gV 31 + p%V 23 + pgV 32 , (7) 

where fi^ v are functions of the Rabi frequencies and re- 
laxation constants and cr^ v are in the absence of in- 
teractions. In practice, one would need to know these in- 
tegrals explicitly. This would require to calculate atomic 
correlators, which is difficult because n-body correlators 
depend on (n + l)-body correlators. 

There are various approaches how to close the system 
of equations for the correlators. One possible way is to 
expand the correlators in powers of the probe Rabi fre- 
quency ftp [23]. Here, we are mostly interested in the 
first two terms in the expansion of <r 12 and cr 33 
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Similarly, all other averages can be factorized in a leading 
term in Q p , and power series in |^ p | 2 . It is easy to 
conclude that Vis °° |^p| 3 an d ^23 ~ |^p| 4 - Therefore, 
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V23 and V32 do not contribute to a\ 2 . We will use the 
propagation equation for (d~l 3 d- J 33 } to illustrate how the 

(3) 

system of equations to calculate a\ 2 gets closed 
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Note that all terms in the sum of three-body averages are 
proportional to \ft p \ 5 in the lowest order except the term 
corresponding to s = j, which is proportional to |^ p | 3 . 
Therefore, only this s = j term contributes to the third- 
order nonlinearities. In a steady-state, the 
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in Eq. ( 10 ) satisfy the following equation 
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We see that this equation does not contain any three- 
body terms. Similarly, the equations for the other cor- 

f3) 

relators ocr^ a p also have no three-body terms and thus 

f3) 

a 12 i s f u lly determined by the expansions of the two- 
body correlators and single atom averages. Even though 
there are many two-body correlators, it turns out that we 
do not need to handle all of them at once because several 
smaller independent subsets of equations are obtained af- 
ter making the expansion. In Appendix [Aj we give more 
details on the calculation of two body correlations using 
the expansion in powers of Q p . 

For stronger excitation, many expansion terms are 
needed in order to get convergent results so this approach 
is not convenient anymore. In this case, it is better to 
know approximately the whole sum than to know ex- 
actly a few expansion terms. Our method presented in 
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Appendix [B] estimates the collisional integrals without 
making use of the expansion. Because of this, it can give 
reasonable, albeit approximate, results in a larger range 
of Q p than a truncated expansion with a few first exact 



terms (some examples are given in Section IV). 

In more details, the calculation proceeds as follows. 
The equations contain various local correlators and non- 
local collisional integrals, that are a consequence of inter- 
actions which couple different atoms and make the calcu- 
lation a difficult many-body problem. In our method we 
first slightly modify the ladder approximation [24] which 
is commonly applied to three-body collisional integrals in 
order to substitute them with two-body ones. The small 
modification is that just before applying the ladder ap- 
proximation we split the collisional integrals involving 
three atoms in two parts corresponding to two different 
leading dependences in power of fl p 



ksi&aS&h&iv) * Sspkji&aS&iv) + (12) 

One may argue that (V a s — hi^ou 0^^3)0 ^ v should be used 
in the last equation instead of V^cr^, but the difference 
between these two expressions is insignificant for small 
033. By conjugating this equation one gets a similar ex- 
pression for the correlators cras a ^ u . The approximation 
(12) is exact up to the lowest order in powers of Vt p in 
winch the interaction dependence starts to appear. For 
instance, it is the third-order for 0013,33 but the fourth 
order for 0023,33 and 01733,33. 

In the second step, we obtain an algebraic system that 
only contains the collisional integrals of interest, as shown 
in Appendix [Bj Since we only need four V^ v integrals to 
determine all a a p, we are ultimately solving four second- 
order polynomials in V^ v . This is certainly much simpler 
than calculating a large number of coupled many-body 



equations. After solving this system for four V^, one 
can get all single-atom averages using the solution ^ of 
Eqs (§-([6]). 



III. SCALING OF BLOCKADED SINGLE- ATOM 
AVERAGES 

In this section we establish some relationships between 
two ways of describing the effects of interactions in Ry- 
dberg ensembles. The first one is more intuitive and it 
is based on the concept of excitation-blockade while the 
other approach is more formal and it uses collisional in- 
tegrals V a p to evaluate the effects of interactions. 

The concept of an excitation blockade sphere is very 
useful to get an intuitive physical picture of the effects 
of Rydberg-Rydberg interactions on <ti 2 . In this pic- 
ture, each Rydberg atom prevents surrounding atoms 
from being excited into Rydberg states and thus these 
rib atoms behave effectively as two-level atoms. Conse- 
quently, a\ 2 of these blockaded atoms is equal to 0i 2 ev of 



two-level atoms. On the other hand, a\ 2 of unblockaded 
atoms is just <jf l 2 ev of noninteracting atoms. Therefore, 
<Ti2 in an ensemble of N atoms of which N r are Rydberg 
atoms, is a weighted average of crf 2 ev for (N — N r nb) un- 
blockaded atoms, and al 2 ev for N r rii) blockaded ones : 
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where p r = N r /N = 033. Essentially the same averaging 
was also used in [25] to describe interaction effects in 
microwave dressed Rydberg atoms. In [25 , the radius of 
the blockade sphere was defined by the condition that the 
Van der Waals shift is equal to the two-photon resonance 
width, and was later used to calculate n b as the number 
of atoms in a blockade sphere. Also, p r was replaced by 
the Rydberg excitation probability £>3 of noninteracting 
atoms, which is acceptable if the blockade spheres do not 
overlap. Another form of universal scaling supported by 
Monte-Carlo simulations is presented in [26] . 

Here, we want to use Eq. ( [13] ) as a definition of n& 
because we can calculate 012 in the presence of interac- 
tions for low probe-light using the rate equations. From 
this perspective, is generally a complex number. Sim- 
ilarly to Eq. (13), it is easy to define a mean value of 
any for /i, v 7^ 3 but extending it to 0^3 and 03 M 
is a bit less straightforward. However, one may expect 
a similar result for 033 = p r « p 3 (l — n b ps) which just 
means that n^p^N atoms are not available for Rydberg 
excitation. Obviously, the n b corresponding to 033 has 
to be real. Therefore, it is also interesting to determine 
to what extent the same parameter can characterize 
the interaction effects on 012 and 033. 

The change in 012 induced by interactions is a^ 11 — 



012 — 0"i2 • To derive in Eq. (13), we start with the 



solutions ([7]) which are now functions of the collisional 
integrals V a p. We take the low probe-light limit of these 
solutions since we know exactly the collisional int egrals 
in this case. Instead of using the exact expression (All) 
for the CRT^a/3 correlator, which is too complicated for 
the theoretical analysis in this section, we use the much 
simpler relation (A12), which is still a very accurate ap- 



proximation in the dispersive regime. 

We assume that the probe field is sufficiently low so 
that the lowest correlation orders correctly describe the 
effects of interactions. Using the solutions Eq. ^ for 012 
and expanding it in powers of Q p , we find the relation 
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where is the lowest expansion coeffi cient of Vis = 
fi p E s= iVi ( 3 2s+1) |fip| 2s . According to Eqs. JXlil, v}^ is 
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where the spatial integral i& depends on the atomic den- 
sity 77, on the Rydberg potential kij and on the effec- 
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tive relaxation constant T « Ti 3 + Vt 2 /Ti2 of the Ry- 
dberg transition, light-shifted and power-broadened by 
the strong control beam: 
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Combining Eqs. (14) and (15) we get 
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Direct inspection of Eqs. (13), (17), and (18) yields the 
expected identification n b - 
Waals interaction k = —Cq/v 



- I b , and for the Van der 
6 we find 



n b = I b = 
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Strictly speaking, n b is a complex number and can even 
be dominant ly imaginary. However, using the atomic 
cloud as as dispersive optical medium puts constraints 
on the signs of A2, A3 and Cq. In order to avoid ab- 
sorption, one must keep all atoms far-detuned from all 
resonances. One consequence of this is that the power- 
boadened linewidth Re(T) « 713 + £7;? /A?, of the Rydberg 
transition must remain smaller than its light-shifted ef- 
fective detuning — Im(T) ~ A3 — Q 2 ,/ for strong con- 
trol beams, this requires A 2 and A 3 to have opposite 
signs. In addition, this effective detuning must have the 
same sign as Cq, otherwise some atomic pairs will be 
resonantly excited into Rydberg states and the optical 
losses will strongly increase. As a consequence, in the 
dispersive regime iT/Cq is dominantly real and positive, 
so that according to the relation (19) n b is dominantly 
real and positive as well: 
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In contrast, if the effective two-photon detuning and Cq 
have opposite signs, n b is dominantly imaginary. Using 
eq. (20) and previous ones, one recovers the formulae 
given in the Appendix of ref. [13]. 

From an experimental point of view, it is also worth 
noticing that in ref. [13] . the control beam is actually a 
standing wave, so that the control Rabi frequency is span- 
ning the whole range from zero to its maximum value. In 
order to keep the above conditions valid in all points in 
space, it is not possible to stay on the two photon res- 
onance as considered in ref. [23], and one has to use a 
small but non-zero value of A3. 



Now we want to find out how 033 is affected by inter- 
actions in the low probe-field limit. However, since 
is the lowest term affected by interactions, all four colli- 



sional integrals contribute to 033 = 033 — 033 . In this 
(16) case we can define the scaling parameter h b as follows 
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This h b is real and has to be positive in the blockade 
regime so clearly n b and n b cannot be exactly the same. 
We can compare these two scaling parameters only in the 
dispersive regime, where n b is also dominantly real and 
positive. Even though one formally needs to calculate the 
forth-order nonlinearities to account for the effects of in- 
teractions on 033^ , this nonlinearity is only due to binary 
interactions, which means that our method in Appendix 
[b] calculates exactly. By expanding Eq. (J^Fj) for 0-33, 
we find that the contributions from V13 and V31 to 033^ 
are approximately 
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Our calculation of V23 for very low probe laser intensity 
shows that the contributions from V23 and V32 are much 
smaller than those from V13 and V31. Therefore, both 
n b and n b depend only on V13 (^31) in the low probe- 
light limit. In this limit we can use the simple relation 
V13 = iQ p Q c psn b /Ti2 in the expression (22) to get the 
following relation 
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One can verify that £1 w 1 but £2 is slightly larger than 
6 for the parameters of Ref. [13] with A3/7 = 1/3. This 
means that the difference between Re[n&] and h b can be 
rather large unless Im[n&] is very small. Hence one can 
write Re[nfr] ~ n b ~ fi b only for large enough two-photon 
detuning, e.g. A3/7 = 2, as will be also apparent from 
the numerical results shown on Fig. [4] 



IV. NUMERICAL RESULTS 

In this section we present our numerical results in the 
dispersive regime. We use physical parameters very sim- 
ilar to the peak values in experiment [13] performed with 
an ensemble of 87 Rb atoms: 77 = 0.04 /im -3 , A2/7 = 
-25, 713/7 = 1/10, and 7/2tt = 3 MHz. The single- 
photon detuning A 2 was chosen to be negative in order 
to avoid heating the atoms with the probe beam, which 
means that dispersive effets could be observed for positive 
two-photon detunings A3 > and attractive Rydberg 
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states Cq > (see section III). Therefore, we focus on 
attractive Rydberg states |3) = \nD 5 / 2 ,mj — 5/2) with 
principal quantum numbers n = 46, 50, 56, and 61 as 
used in Ref . [T5] , where the numerical values of Cq/j are 
respectively 2400, 5000, 15000, and 36000 fim 6 . Assum- 
ing that the control laser irradiance remains constant, we 
take f2 c /7 = 3 for the state n = 50 and use the n -3 / 2 
dependence of the transition dipole to calculate Q c for 
the other states. 

We make two additional assumptions which simplify 
the calculations compared to the actual situation in 
Ref. [13] and still provide a good qualitative agreement 
with the experimental data. First, we disregard the (lon- 
gitudinal and transverse) spatial and temporal variations 
of the Rabi frequencies and of the atomic density, and 
consider them as static and uniform: taking these varia- 
tions into account is not a trivial problem and the quan- 
titative modeling of the experimental results will be pre- 
sented elsewhere. Second, we assume like we did be- 
fore that the interactions between Rydberg atoms are 
described by a single potential kij , which is actually only 
true for S Rydberg states. In our case the nD 5 / 2 states 
with the maximal projection rrij = 5/2 were targeted. 
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This implies that all ungerade nD$/ 2 
curves of Rb should be considered [27 j. All these poten- 
tials are attractive but very different in strength, ranging 
over two orders of magnitude. Since we are interested in 
the dispersive regime where none of these potentials be- 
comes resonant, and since we consider only global prop- 
erties obtained after radial and angular integration over 
the atomic cloud, we can try to reproduce their overall 
effect with a single effective potential [28] |29] . How this 
potential is defined in detail depends on the problem and 
on the purpose for which it is used, and our procedure 
will be explained in detail in a separate publication. It 
turns out that the values for our effective Cg ff differ by 
less than 20% from those in [28 . The difference arises 
because we average \/C§ over the relevant potentials to 
get an effective \J Cg ff , which is the quantity of interest 
in our problem, while in [28] Cq was averaged to get Cg ff . 

The optical response of the cloud to the probe beam 
is characterized by its optical susceptibility X12, propor- 
tional to cr 12 /O p . By using appropriate normalization 
this proportionality constant can be disregarded, and 
in the following we will call "susceptibility" the ratio 
<7i2/^ p . We can define the following normalized suscep- 
tibility S norm as a measure of nonlinear effects induced 
by the Rydberg interactions 



FIG. 2. Normalized susceptibility ([26} for n = 46, 50, 56, 
and 61. For noninteracting systems S = 1 holds everywhere 
while S = is reached in the fully blockaded regime. The 
slopes (dashed lines) at the origin are given by the third- 
order interaction-induced nonlinearities and they follow the 



n ^5/2 potential \[C§ oc n 11 ^ 2 scaling as a function of the principal quantum 



number n. The solid curves are calculated using the method 
in Appendix [B] 



shifts $ of the cavity resonance frequencies, so it is useful 
to consider the normalized real part of the susceptibility 
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because, in a single- mode cavity, it is equal to the nor- 
malized cavity resonance shift 4>(|fi p | 2 )/3>(0) which has 
been measured [15] . 

On Fig. [2] the full curves present the shift S calculated 
using the methods described in Appendix [Bj whereas 
their initial slopes (dashed curves) are determined by the 
third-ord er no nlinearities evaluated using the exact ex- 



pression (All) for crcr 13 33 , 
because 
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In the absence of interactions, (j\ 2 = crf 2 ev holds so that 
1. In a fully blockaded system g\ 2 ~ cr1 2 ev so 



u norm 

that S, 



vanishes. If the atoms are placed in an optical 
cavity, all the coherences in Eq. ( [25] ) should be substi- 
tuted by their projections onto the cavity mode(s). In 
this case the nonlinear phase shifts are mapped into the 



As pointed out in Appendix [B] our approximate method 
should reproduce the third-order nonlinearities exactly 
and this is manifestly true in Fig. [2] Since the actual 
normalized susceptibilities rapidly deviate from the ini- 
tial third-order dependence, all the more as n is large, this 
figure also shows that the higher order (> 3) nonlineari- 
ties are important and, in fact, quickly become dominant 
in our parameter range. Even though all the curves show 
signs of saturation, they have not yet reached the fully 
blockaded regime (5 = 0), which is also confirmed by the 
experimental results [13]: the normalized susceptibilities 
are getting saturated before reaching the fully blockaded 
regime and then slowly approach the limit 5 = 0. 
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FIG. 3. Real and imaginary part of susceptibility G\il^t v 
in the vicinity of the two-photon resonance. All parameters 
are the same as in Fig. [2] for n = 61, except that Q p is 
fixed: \Q P \ 2 — 7 2 /2. The susceptibilities in the noninteracting 
(blue) and interacting case (red) are shown together with the 
susceptibility of two- level atoms (dashed green). In the inset 
plots, the noninteracting erf l 2 ev /il p (blue) is compared with the 



truncated expansion + |^p| 2 o"i2 j (dashed black) 



2 (3) 



In Fig. [3] we show the real and imaginary part of 
<7i2/Hp in the vicinity of the two photon resonance for 
n = 61. The physical parameters are the same as in 
Fig. [2] for n = 61, except that we take the maximal 
l^pl = t/v^2 considered in Fig. [2] To compare these 
results with the two limits of ^norm, we show cr 12 /O p 
for non-interacting and interacting atoms together with 
G\2/&p of two-level atoms. Even though the strongest 
interaction effects in Fig. [2] are seen for n = 61, the 
corresponding cr 12 has not yet reached the fully block- 
aded limit (Ji2 V ' As A3 approaches zero, the interaction 
and blockade effects become more pronounced. In the 
inset plots, we show the noninteracting af l 2 ev /ft p and its 



truncated expansion -\-\Q p \' 2 a[°2 corresponding to the 
third-order approximation of (T12. We see that this third- 
order approximation completely fails for negative detun- 



2^(3) 



FIG. 4. Scaling parameters Re[nb] (solid) and hb (dashed) for 
A3/7 = 1/3 (red), A3/7 = 1 (green), and A3/7 = 2 (blue) as 
functions of probe-laser intensity. All the other parameters 
are the same as in Fig. [2] for n — 50. 



ings even for noninteracting atoms. For A3/7 = 1/3, 
which was used in Fig. [2j it remains satisfactory in the 
noninteracting case, but in the presence of interactions 
the range of Q p where this expansion remains valid is 
dramatically reduced. This confirms the result of Fig. [2j 
where the actual curves quickly deviate from their initial 
slopes given by the third-order non-linearity. The block- 
ade spheres are expected to be "stiff' but not really hard 
balls. Therefore, it is expected that n b slowly decreases 
as Q p increases. To take this into account, we can use 



Eqs. (13) and (21) as definitions of and for any Q p . 
Since is complex and h^ is positive, we can only com- 
pare them if they are both positive. In Fig. |4j we show 
how Re [ni,] an d depend on the probe-laser intensity for 
A3/7 = 1/3, 1, and 2. All the other parameters are the 
same as in Fig. [2] for n = 50. The results demonstrate 
that as A3 increases the imaginary part of 715 vanishes, 
and both scaling parameters and their difference decrease 
as expected. One can also see a significant difference be- 
tween Re[nfr] and for A 3 = 7/3, which is due to the 
rather large factor £2 = 6.4, and to the non-negligible 
(although small) Im[n^]. As for the contributions of V23 
and V32 to we verified that they remain below ~ 3%. 



V. CONCLUSIONS 

In this work we investigated strong dispersive optical 
nonlinearities in ensembles of Rydberg atoms under EIT 
conditions. Using the expansion in powers of Q p , we eval- 
uated the steady-state third-order susceptibility, domi- 
nant for sufficiently low probe fields. In this regime we 
established precise relationships between an intuitive de- 
scription based on Rydberg blockade and a more formal 
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one based on the rate equations and collisional integrals. 
We introduced scaling parameters that encapsulate 
the effects of interactions on single-atom average popu- 
lations and coherences. Generally speaking, the that 
correspond to coherences are complex and can even be- 
come dominantly imaginary for resonant excitation of 
atomic pairs. However, when the Rydberg excitation 
blockade is used in the dispersive, off-resonant regime, 
we found that these scaling parameters become positive 
and very similar to each other, taking the usual interpre- 
tation of the number of atoms in a blockade sphere. 

Beyond the lowest-order regime, like in any many-body 
calculation, the complexity arises from the treatment of 
nonlocal interaction-dependent quantities. We resolved 
this problem by approximating many-body correlations 
by two-body ones, and finding a closed set of algebraic 
equations satisfied by the latter. This method reproduces 
exactly the lowest-order interaction-induced nonlineari- 
ties, allowing us to determine the range of probe laser 
intensities where the lowest-order expansion is valid. As 
expected, this range decreases with increasing Rydberg 
interactions. We calculated the optical response of the 
atomic cloud under conditions similar to those in a re- 
cent experiment [13], and we justified the theoretical ap- 
proach that was used to interpret the results of the exper- 
iment. The numerical results are in very good qualitative 
agreement with the experimental data. A full quantita- 
tive analysis requires to take into account the spatial and 
temporal variations of Rabi frequencies; this will be the 
object of another work. 
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l(T 12 + r 23 )<TCri2,23 = ^23 ~ ^13,12 + 



(A3) 
,(A4) 



•m i \ (3) (2) (2) 

l(Tl3 + 722)^13,22 = ^13,21 - ^13,12 



• /-p* I "j— * \ (3) (2) . (2) 

H r 23 + r i2)^12,32 = ^31,12 + a 32 



n c [o-o-{ 3 2 ] 22 - o-a{ 3 2 ] 33 + vo-fl 23 ] , ( A5) 

(2) (2) 
= ^13,21 - ^13,12 + 

[°v{ 3 2 ] 22 + <*7^23 - a(J i%2] > ( A6 ) 

(A7) 
(A8) 



7 (3) 
7 12,33 



.(3) 
7 13,32 



7 12^22] 



Vn , N (3) (2) (2) . (2) 

H r i2 + 722)^12,22 = acr 12,21 - ^12,12 + a 22 



n c [aa[ 3) 



2,23 



7 (3) 
7 13,22 



7 (3) 1 
7 12,32j ' 



All expansion terms are solved order by order so it is as- 
sumed that the first- and second-order expansion terms 
have been solved before solving the third-order terms. 
Therefore, all second order averages in Eqs. (A1)-(A8) 



are considered to be known functions. Only two corre- 



lators cra^ 333 and o-a\° 3 23 in the first two equations of 



.(3) 



(3) 

kijO~o~ g0 uu . We can eliminate all the other correlators in 



the system (A1)-(A8) are related to interaction terms 



Eqs. (A1)-(A8) to get 



(F 1 -k ij )aa {3 3 ] 33 - K x gg^ 23 

.(3) _ j^O) _ 



(3) _ Aikj. 



P H~ kijQ 

(F 2 -k ij )ou\t> 23 - K 2 ou£> 33 = ^ + G 2 , (A10) 



P H~ kijQ 



where all F a , K ai G ai P and Q are lengthy but known 
functions of control Rabi frequency tt c and relaxation 
constants T^, 7^. From eqs. ( A9|A1Q ), we can find a 



general expression for the radial dependence of acr 



(3) 

13,33' 



Appendix A: Third-order solution 



Here we show the main steps which yield the exact so- 
lution for the third-order susceptibilities in the presence 
of interactions. Because this solution is too long to be ex- 
plicitly written, we will just provide its general form. It 
turns out that the rather large system for two-body cor- 
relators decouples into several independent subsets which 
can be solved order by order. We are only interested in 
the correlators that are relevant for the collisional integral 
V13. Therefore, we end up solving the following system 
(for 733 = 0) 



(ir 13 

(^13 



&ij)ov 13 33 - 



_(3) 
7 12,33 



r (3) 
7 13,23 



.(3) 



'13,32 



,(A1) 



iT 



23 



/r-Wrr (3) 

^lj)OT7 13 o 3 



r( 2 ) 



13,13 

[^12^23 + 0(7 13 22 ~ 0,7 13^ 33] > (A2) particles so that era, 



aa 



(3) 
13,33 



(All) 



/ (r,^)+A(r,^)^+/ 2 (r,^)^ 
fto(r,n)+hi(r,n)fey+h 2 (r,n)fc?.+h 3 (r,n)fc?/ 

where fi and hi are functions of relaxation constants and 
Rabi frequencies. This general form can be analytically 

(3) 

integrated to get the collisional integral V-± 3 which is im- 
portant for the interaction induced optical nonlinearities 



and au 



(3) 

12,33* 



We use the exact based on Eq. (All) when we 



want to precisely know its numerical values, such as the 
slopes in Fig. [2] However, for some theoretical analy- 
sis, it may be more convenient to use approximate but 
simpler expressions. Since we mainly consider the disper- 
sive regime here, we can assume that IA2I is much larger 
than any other frequency. For very large interatomic 
separations Rij there should be no correlations between 
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Rij — » oo. For finite separations R^ and |Ti2 1 ^ the 
correlator is 



T (3) 



(1)^(2) 



0013,33 ~ a 13 a 33 y + i& 



(A12) 



where T w Ti3 + l^/r^ is the effective damping rate 
of the Rydberg transition. This is in fact an excellent 
approximation for |Ti2| ^> 



and V32. Therefore, i? m and 3? n are at most second or- 
der polynomials of the four collisional integrals. To ad- 
ditionally clarify those P(Q)-correlators we give a few 
examples. For instance, <t<ji3 5 33 is a P-correlator while 
0022,33 is a Q-correlator whose corresponding 3? n is equal 
to zero. On the other hand, au-22,13 is a Q-correlator with 
nonzero 5R n . 

Using Eqs. (B3), we can eliminate all Q s from Eqs. 



(B2) as follows 



Appendix B: Collisional integrals 

In order to evaluate the single-atom averages in the 
presence of interactions, we need to know the collisional 
integrals V^ 3 = Ej^i kij (0^3033)- There are several well 
known difficulties of doing this. The first one is that 
(n + l)-body collisional integrals are needed to evaluate 
n-body correlators. We use a slightly modified ladder 
approximation (12). The ladder approximation is a stan- 



dard approach to this problem. Another problem is that 
the steady-state equations mix local correlators, as func- 
tions of interatomic separations between involved atoms, 
and nonlocal collisional integrals. Our goal in this ap- 
pendix is to explain how one can get a closed algebraic 
system that contains only two-body collisional integrals, 
which is easy to solve numerically. 

We start with the rate equations for two-body corre- 
lators (J(Ja(3,nv assuming that the approximation (12) has 
been utilized. Most of these equations do not explic- 
itly contain the terms kijoa^^, where k^ = —C s /Rfy 
(From now on we use k and R instead of k^ and Rij.) By 
P m , we denote two-body correlators (JOafi^v whose time 
derivative contain the term k P m and Q n are all other 
two-body correlators. As mentioned previously, we as- 
sume the symmetry ou a p^v = o&^ap. In addition, all 
substitutions 



0011, fiu = 0/izy — 0022,/i^ — 0033,/ii/ 



(Bl) 



are made. This results in ten equations for ten P m cor- 
relators and twenty six equations for twenty six Q n ones 
which are conveniently written as follows 

kPm ~ ^ ] &mqPq H~ ^ ] ^msQ s H~ Rm-> (B2) 

q s 



where the coefficients a mq , 6 ms , c n 



d nq are either 



the Rabi frequencies or relaxation constants and i? m 
and 3? n contain terms which are essentially constant 
on the length scale of interaction induced correlations. 
These R m and 5R n are in general linear com binations of 



(originating from the substitution (Bl)) and terms 
0/ii/^a3 and cr^Vsa (originating from the ladder approx- 
imation (12)). According to Eq. (|7|), all are linear 
combinations of four collisional integrals V13, V31, V23, 



d nq )P q + (B4) 

q n 

^2(b ms + ^2 a mnCns)Qs + Rm + ^ Ol mr J& n . (B5) 



Imposing the condition 

n 

for each m and s, we get 

kPm — ^ ^ ttmqPq H~ Rn 



(B6) 



(B7) 



where cL 



Iraq H~ ^n^mn 

d nq . We can introduce ma- 



trices a, 6, c, a as matrices corresponding to the coef- 
ficients a mq , b ms , c ns , d nq , and a mn , respectively. From 
Eq. (B6) we obtained a = — 6c _1 . The matrix form of 
Eq. fcj) is 



fcP = (a-&c _1 d)P + R, 



(B8) 



For any eigenvector A and corresp ondi ng eigenvalue A of 
the matrix (a — &c _1 <i) T , relation (B8) implies 



fcA.P = AA.P + A.R, 



and consequently 



A.P 



A.R 

k- A 



(B9) 



(BIO) 



Note that A.B is not a sca lar product but just the sum 
S m A m P m . From Eq. (BIO), we see that each eigenvector 



A imposes one condition on the collisional integrals 
V = 7?y 4 d 3 RkP 



as follows 



A.V 



r] / d 3 R 



k-X 



\.R = F(\)\.R. (Bll) 



Introducing matrices 





' Ai " 




' F(Ai)Ai " 


u = 




, F = 










F(Aio)Aio 



(B12) 
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we get 

V = (7 _1 FR. 
For the Van der Waals interaction k = —Cq/R 6 we 



find F(X) = (2tt 2 7 ? /3)V / C67A. 

We do not need to consider all ten collisional integrals 
V but only four V13, V31, ^23, and V32 since any compo- 
nent of R only depends on these four collisional integrals. 
Once these integrals are determined, the single-atom av- 
erages are obtained using the solution ^ of Eqs. ([2])-([6|. 
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